{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/Qiskit/qiskit-tutorials/master/images/qiskit-heading.png\" alt=\"Note: In order for images to show up in this jupyter notebook you need to select File => Trusted Notebook\" width=\"500 px\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _*Measurement Error and Mitigation*_ \n",
    "\n",
    "\n",
    "***\n",
    "### Contributors\n",
    "David McKay and Yael Ben-Haim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "The last step of a typical quantum experiment is to perform a measurement on the qubits in the circuit. Although the qubit state $|\\Psi\\rangle$ (or more generally the density matrix $\\rho$) is the general description of the quantum state, in a typical strong projective measurement our measurement _projects_ the general state into a specific computational state $|x\\rangle$ (where $x$ is a bitstring, e.g.,  `1001010`) The probability of measuring bitstring $x$ is given by:\n",
    "$$P_x = \\mathrm{Trace}(\\langle x|\\rho|x \\rangle)$$\n",
    "Therefore, the measurement process is stochastic. The above distribution of $x$ given a state $\\rho$ is true only in the absence of measurement errors. There are multiple sources of possible measurement error, all of which are dependent on the physical mechanism of measurement in the system. For superconducting qubits coupled to readout cavities [[1](#ref1),[2](#ref2),[3](#ref3),[4](#ref4),[5](#ref5)] the state of the qubit is determined by measurement the response of a microwave tone incident on the readout cavity. The cavity signal is measured for some time where $V(t)$ is the complex amplitude of the signal which is converted to a single complex number based on a measurement kernel \n",
    "$$V = \\int_0^{T} V(t) K(t) dt $$\n",
    "which is then turned into a _bit_ by a nonlinear discriminator [[6](#ref6)]. The simplest example being if $|V|<V_0$ then the qubit was in state 0 and otherwise the qubit was in state 1. \n",
    "\n",
    "As discussed in [[6](#ref6)] there are classical sources of noise on the signal that lead to misidentification of the qubit state, but it can also happen that the qubit decays due to $T_1$ during the measurement. There are other sources of crosstalk (to numerous to enumerate) such as classical crosstalk on the lines and crosstalk between resonantors on chip. All of these issues lead to a new probability distribution $\\tilde{P}_{\\rho}$ for a given state. Given certain assumptions about these errors and appropriate calibration we can attempt to correct the skew in the probability distribution _on average_. \n",
    "\n",
    "\n",
    "### References\n",
    "\n",
    "[1]<a id=\"ref1\"></a> Alexandre Blais, Ren-Shou Huang, Andreas Wallraff, S. M. Girvin, and R. J. Schoelkopf, Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation, https://arxiv.org/abs/cond-mat/0402216\n",
    "\n",
    "[2]<a id=\"ref2\"></a> Jay Gambetta, Alexandre Blais, D. I. Schuster, A. Wallraff, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf. Qubit-photon interactions in a cavity: Measurement induced dephasing and number splitting\n",
    "https://arxiv.org/abs/cond-mat/0602322 \n",
    "\n",
    "[3]<a id=\"ref3\"></a> Alexandre Blais, Jay Gambetta, A. Wallraff, D. I. Schuster, S. M. Girvin, M. H. Devoret, and R. J. Schoelkopf. Quantum information processing with circuit quantum electrodynamics. https://arxiv.org/abs/cond-mat/0612038\n",
    "\n",
    "[4]<a id=\"ref4\"></a> Jay Gambetta, W. A. Braff, A. Wallraff, S. M. Girvin, R. J. Schoelkopf. Protocols for optimal readout of qubits using a continuous quantum nondemolition measurement. https://arxiv.org/abs/cond-mat/0701078\n",
    "\n",
    "[5]<a id=\"ref5\"></a> Jay Gambetta, Alexandre Blais, M. Boissonneault, A. A. Houck, D. I. Schuster and S. M. Girvin. Quantum trajectory approach to circuit QED: Quantum jumps and the Zeno effect. https://arxiv.org/abs/0709.4264\n",
    "\n",
    "[6]<a id=\"ref6\"></a> Colm A. Ryan, Blake R. Johnson, Jay M. Gambetta, Jerry M. Chow, Marcus P. da Silva, Oliver E. Dial and Thomas A. Ohki. Tomography via Correlation of Noisy Measurement Records. https://arxiv.org/abs/1310.6448"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code imports\n",
    "=============="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-18T15:13:54.475305Z",
     "start_time": "2018-12-18T15:13:53.026353Z"
    }
   },
   "outputs": [],
   "source": [
    "# Import general libraries (needed for functions)\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "# Import Qiskit classes\n",
    "import qiskit \n",
    "from qiskit import QuantumRegister, QuantumCircuit, ClassicalRegister, Aer\n",
    "from qiskit.providers.aer import noise\n",
    "from qiskit.tools.visualization import plot_histogram\n",
    "\n",
    "# Import measurement calibration functions\n",
    "from qiskit.ignis.mitigation.measurement import (complete_meas_cal, tensored_meas_cal,\n",
    "                                                 CompleteMeasFitter, TensoredMeasFitter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constructing a Full Calibration Matrix\n",
    "\n",
    "The assumption of the error mitigation technique is that we can prepare each of the basis states with very low error. Given this assumption, in separate experiments we can prepare one of the $2^n$ states and then measure the outputs in all $2^n$ states ($n$ denotes the number of qubits). Normalizing these outputs and making each set of output probabilities for a given prepared state the columns of a matrix we obtain the matrix $\\mathbf{A}$ which translates the ideal probability distribution of the state $\\rho$ ($P_\\rho$) into the experimental probability distribution $\\tilde{P}_{\\rho}$\n",
    "$$\\tilde{P}_{\\rho} = \\mathbf{A} \\cdot P_{\\rho}$$\n",
    "\n",
    "\n",
    "**Code**\n",
    "\n",
    "The code below constructs the calibration matrix for 2 qubits with error artificially put into the Aer simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ideal calibration matrix:\n",
      "[[1. 0. 0. 0.]\n",
      " [0. 1. 0. 0.]\n",
      " [0. 0. 1. 0.]\n",
      " [0. 0. 0. 1.]]\n",
      "Noisy calibration matrix:\n",
      "[[0.81  0.244 0.201 0.045]\n",
      " [0.087 0.646 0.032 0.192]\n",
      " [0.092 0.023 0.693 0.213]\n",
      " [0.011 0.087 0.074 0.55 ]]\n"
     ]
    }
   ],
   "source": [
    "# Generate the calibration circuits\n",
    "qr = qiskit.QuantumRegister(2)\n",
    "meas_calibs, state_labels = complete_meas_cal(qubit_list=[0,1], qr=qr)\n",
    "# Generate a noise model for the 2 qubits\n",
    "noise_model = noise.NoiseModel()\n",
    "for qi in range(2):\n",
    "    read_err = noise.errors.readout_error.ReadoutError([[0.9, 0.1],[0.25,0.75]])\n",
    "    noise_model.add_readout_error(read_err, [qi])\n",
    "backend = qiskit.Aer.get_backend('qasm_simulator')\n",
    "job_no_noise = qiskit.execute(meas_calibs, backend=backend, shots=1000)\n",
    "job_w_noise = qiskit.execute(meas_calibs, backend=backend, noise_model=noise_model, shots=1000)\n",
    "cal_results = job_no_noise.result()\n",
    "meas_fitter = CompleteMeasFitter(cal_results, state_labels)\n",
    "print(\"Ideal calibration matrix:\")\n",
    "print(meas_fitter.cal_matrix)\n",
    "cal_results = job_w_noise.result()\n",
    "meas_fitter = CompleteMeasFitter(cal_results, state_labels)\n",
    "print(\"Noisy calibration matrix:\")\n",
    "print(meas_fitter.cal_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that with noise when we prepare the state $|11\\rangle$ and measure we get counts in states other than $|11\\rangle$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAE6CAYAAABeVIXiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxV1bn/8c+ThFlAgmlIghAGDWEwBrRXqeKAKKhVC9jaWrVaoWrFAf2p1VbxtmCl2mrVXi3e1ulSh4J1pHIVr4hDbSNGhhjBSBAIGAMCMiUkz++PfZKGkIQcODk7Jt/363VeOWfttdd5NkvMw1prr23ujoiIiIjEV0LYAYiIiIi0RUrCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKQFHYAYerZs6f36dMn7DBERKQF+eCDD75w95Tqz3l5ed9ISkp6GBiKBi+k6aqApbt37750xIgRn9dXoU0nYX369GHBggVhhyEiIi1IcnJyce3PSUlJD/fq1Ss7JSVlU0JCgjbXlCapqqqy0tLSwevXr38YOKu+OsroRUREGjc0JSVlixIwiUZCQoKnpKRsJhhBrb9OHOMRERH5OkpQAib7I/LfTYO5lpIwERERkRAoCRMREREJgZIwERERkRAoCRMREZFmMXDgwCEvvvhi15bU7tSpU9OnTp2aHuuY9keb3qJCRERkf0y6hxHN2f6sa8hrat2MjIxhZWVl7RISErxTp05VJ5544uY//elPq7t3717VnDE2xcqVK5e19HZnzJiRMnv27EM+/vjjTt/+9rc3zpkzZ1VTj2dkZAxbuHDhRwMGDKjYn+/WSJiIiMjX3JNPPrli+/bti//1r38t//DDD7v87Gc/S6tbp6Jiv/KE/RLP7zpQGRkZFTfeeGPJueee+0W0x8eMGfPlM888c/D+freSMBERkVaiX79+FSeffPLmgoKCThCM1Nxyyy29Dj/88MFdunQZXlFRwapVq9qddtppA3r06JGTkZEx7Fe/+tU3qs/PyMgY9rOf/azXgAEDhnTr1u3IiRMnZm7fvt2qj9988829Dj300KFdunTJHTBgwJDHHnvs4Nrn1v2ujIyMYX/729+6Vh+/9dZbUw8//PDBXbt2PfKMM87oX7vtRYsWdc7Ozh7cpUuX3HHjxvU/44wz+l911VX1ThvWbrcpbTfmoosu+vKCCy74smfPnrujPX7OOed8+cILLygJExERaetWrlzZ7rXXXut+xBFHbK8umzNnTvLLL7+8YuPGjYsTEhI444wzBg4bNmx7SUnJh//7v/9b+OCDD6bOmTOnW3X9v/71rz1feeWVj1esWLGkqKio40033VQzqjZw4MBdb775ZuGWLVsW33TTTet+8pOf9CsuLm5X33e1a9eOup599tnk+fPnr1i5cuWSgoKCTvfff/8hADt37rTvfve7A37wgx98sXHjxg/OO++8jfPnz48quWmo7eY0bty4rwoKCjqXlZUl7s/5SsJERES+5n7wgx8M7Nq165GjRo0adMwxx2ydPn16SfWxyy67bMPAgQMrDjroIH/jjTe6bNy4Memuu+4q6dixow8ePLj8ggsuKP3LX/6SXF1/0qRJnw8cOLAiNTW18sYbbyx59tlna45dcsklmzIzMysSExOZNGnSpr59++568803u9T3XfXFefnll2/IzMysSE1NrTz11FM3f/DBB50AXn/99S67d++2W2655fMOHTr4RRdd9OURRxyxLZo/g4babk4dOnTwb33rW1tqJ7HR0MJ8ERGRr7nZs2evPOecc7bWd6xv3741C7SKioral5aWtu/ateuR1WVVVVV21FFH1Zzbp0+f8ur3AwYM2FVaWtq++vP999/f8/77709du3Zte4AdO3YklpaW1uQStb+rPunp6TXHO3fuXFVSUtIO4LPPPmuXmppakZCQULtueT1NRN12cxs6dOiOJUuWdAY2RXuukjAREZFWzMxqRqUyMzPLMzIydhUXFy9tqP7q1atrkq6ioqL2KSkp5QAff/xx+6lTp/Z94YUXPh49evRXSUlJDBo0aLD7vwe9an9XNDIyMio2bNjQrqqqiupEbN26de379eu3a3/ai6cFCxZ0u/XWW9ftz7majhQREWkjTjzxxG1dunSpvOWWW3p99dVXtnv3bv75z392fOONNzpX13n44YdTPvnkk3YbNmxIvPPOO9POPvvsTQBbt25NMDN69epVAXDvvff2XLlyZUym/EaPHr0tMTHR77jjjm9UVFTwxBNPHPzhhx922feZB66iooLt27dbZWWlVVZW2vbt26323Z2NHS8tLU385JNPOp122mlf7c93ayRMREQkStHs49WSJCUl8fLLL6+cMmVK78zMzCPKy8utX79+O2+//fa11XUmTJiw8dRTTz38888/bzdmzJgv77jjjhKAESNG7Jw8efKGUaNGZSckJPiECRPKcnNz9yv5qKtjx47+1FNPfTJ58uTM6dOnZ5xwwgmbTzrppM0dOnRo9gen33jjjem/+93vam4+6NKlS/K1115b8tvf/nbdvo4/88wz3UeNGrU5KWn/0imrPYzY1uTm5vqCBQvCDkNERFqQ5OTkPHc/qvpzfn7+qpycnHr3kGptMjIyhj3wwAOrGlpfFk9HHHHEoB//+MelV199dVks263eLb86yToQ48aN6//9739/44UXXvhlQ3Xy8/MPycnJyazvmKYjRUREJHQvvfTSQatXr06qqKjgvvvu6/nxxx93Puecc7aEHVdjOnXqVHUgMWo6UkREREJXUFDQ8cILLxywY8eOhN69e+965JFHPtnX3Zb74+STT47ZKN/cuXNXHcj5cU/CzOwK4P8BacAy4Bp3f7OR+u2BnwMXAOnABuAud/99rToTgF8CA4BPgFvc/dlmuwgREZFWaO3atUvC+u7rr7/+i+uvv77Zp33PPPPM0Kdaq8V1OtLMvgfcC8wAcoG3gXlm1qeR054ExgKTgSzgXODDWm0eCzwF/A9wZOTnM2b2H81xDSIiIiKxEO+RsKnAI+4+K/J5ipmNBS4Hfla3spmdCowGBrh7dXa8qk61a4DX3X165PN0MzspUv79GMcvIiIiEhNxGwmLTCuOAObXOTQfGNnAaecA/wSmmtkaM1thZr83s4Nq1Tm2njZfaaRNERERkdDFcyTsECCRYE1XbRuAUxo4pz9wHLALmAAcDNxHsDZsYqROrwba7FVfg2Y2mWBqk7S0NN5//30A0tPT6dy5MytXrgSge/fu9O/fn8WLFwOQmJhITk4OhYWFbNsWPM4qOzubjRs3smFD8PW9e/emffv2FBUVAdCjRw/69OlDfn4+AO3atWPYsGEUFBSwY8cOAAYPHkxpaSmlpaUA9O3bFzNj1apVAPTs2ZO0tDSWLg02N+7QoQNDhgxh2bJl7NoVbCQ8dOhQSkpKKCsL7uLNzMzE3SkuLgYgJSWFlJQUli9fDkCnTp3Izs5myZIlVG84l5OTw+rVq9m0KXjqQv/+/SkvL2fNmjUApKamkpycTEFBAQBdunQhKyuL/Px8KisrAcjNzaWoqIjNmzcDMHDgQLZv3866dcFdwGlpaXTr1o3CwkIAunbtymGHHcbixYtxd8yM3NxcVqxYwdatwZR9VlYWW7ZsoaSkRP2kflI/tZB+ev7555k+fTpVVVWcd955XHnllXv001tvvcUvfvELDjnkEMyMK6+8kqOPPpodO3Zw6qmnMmjQICorK0lOTuaXv/wlffv2Zc2aNVx++eVs2bKFoUOH8vDDD/Pxxx+H0k/1qKqqqrKEhIS2u6eT7JeqqioDqho6Hrd9wswsHVgLnODuC2uV3wqc7+5Z9ZwzHzge6OXumyNlpxKMdPVy9w1mVg5c6u6P1TrvQmCWu3doLCbtEyYiEp3KykqOPvpo5s6dS3p6OqNHj2bWrFkMGjSops7s2bP54IMPmDlz5l7nH3rooXz22Wd7lV988cWceeaZTJgwgalTpzJ06FAuueSSZr2WhtSzT9jzvXr1GpySkrJZiZg0VVVVlZWWlnZfv3798pycnLPqqxPPkbAvgEogtU55KrC+gXNKgLXVCVhEQeRnH4IRr/VRtikiIvspLy+Pfv36kZmZCcD48eOZN2/eHklYtNydN998k1mzguXC5513HnfeeWdoSVhdu3fvvnT9+vUPr1+/fijaX1OargpYunv37ksbqhC3JMzdy80sDxgDPFPr0BhgTgOnvQWca2YHuXv1oxEOj/wsjvx8J9LGb+q0+XZMAhcRkRolJSVkZGTUfE5PTycvb+8n+Lzwwgu8/fbbDBgwgOnTp9O7d28Adu7cycknn0xiYiLXXHMNZ5xxBhs3bqR79+5UP/olPT29Zsq0JRgxYsTnQL0jGSIHIt53R/4WeNzM3iNIsC4jWN/1IICZPQbg7hdG6s8GfgH82cymEawJuxf4q7t/HqlzL7DQzG4C/gZ8BziJYC2ZiIjE2dixY5kwYQIdOnTgkUce4ac//SnPPfccAPn5+aSnp7Nq1SrOPvtsBg8eTLdu3UKOWCQccR1WdfenCLaO+DnwAUGidLq7V49q9Ym8qut/RbBovzvBXZJPA28Al9Sq8zZwHvAjgv3DLgS+5+7/aObLERFpc9LS0li7tuZZz6xbt460tLQ96iQnJ9OhQ7Ak94ILLuCDDz6oOZaeng4EC+mPO+44PvzwQ5KTk9m8eTO7d+9usE2R1ijuc9vu/gd3z3T3Du4+ovYifXc/0d1PrFO/0N1PdffO7p7h7j9196116vzV3Qe5e3t3z3b3uXG6HBGRNmX48OEUFRVRXFxMeXk5c+fOZezYsXvUWb/+30ty582bx+GHB6tIvvzyy5q7G8vKyvjHP/5BVlYWZsZxxx1XM1r25JNPcvrpp8fpikTCo2dHiohIkyUlJTFz5kwmTpxIZWUl559/PtnZ2cyYMYPc3FzGjRvHH//4R+bNm0dSUhI9evTggQceAKCwsJCpU6eSkJBAVVUVV199dc2C/mnTpnHppZcyY8YMhg0bxg9/+MMwL1MkLuK2RUVLpC0qRESkrrpbVIg0F91qKyIiIhICJWEiIiIiIVASJiIiIhICJWEiIiIiIVASJiIiIhICJWEiIiIiIVASJiIiIhICJWEiIiIiIVASJiIiIhICJWEiIiIiIVASJiIiIhICJWEiIiIiIUgKOwAREQnPDY/2CDuEmJl50aawQxCJikbCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBErCREREREKgJExEREQkBHFPwszsCjP71Mx2mlmemR3fxPOOM7PdZra0TvmPzMzreXVsnisQEREROXBxTcLM7HvAvcAMIBd4G5hnZn32cV4P4DHgtQaqbAfSar/cfWes4hYRERGJtXiPhE0FHnH3We5e4O5TgBLg8n2c99/Ao8A7DRx3d19f+xXDmEVERERiLm5JmJm1B0YA8+scmg+MbOS8K4BU4FeNNN/JzIrNbI2ZvWhmuQccsIiIiEgzSorjdx0CJAIb6pRvAE6p7wQzGwbcBhzj7pVmVl+1QuASIB/oClwNvGVmOe6+op42JwOTAdLS0nj//fcBSE9Pp3PnzqxcuRKA7t27079/fxYvXgxAYmIiOTk5FBYWsm3bNgCys7PZuHEjGzYEl9S7d2/at29PUVERAD169KBPnz7k5+cD0K5dO4YNG0ZBQQE7duwAYPDgwZSWllJaWgpA3759MTNWrVoFQM+ePUlLS2Pp0mApXIcOHRgyZAjLli1j165dAAwdOpSSkhLKysoAyMzMxN0pLi4GICUlhZSUFJYvXw5Ap06dyM7OZsmSJVRUVACQk5PD6tWr2bRpEwD9+/envLycNWvWAJCamkpycjIFBQUAdOnShaysLPLz86msrAQgNzeXoqIiNm/eDMDAgQPZvn0769ato/rPu1u3bhQWFgLQtWtXDjvsMBYvXoy7Y2bk5uayYsUKtm7dCkBWVhZbtmyhpKRE/aR+Uj81Qz9BD1qLsrKymPSTSLyYu8fni8zSgbXACe6+sFb5rcD57p5Vp34HYDFwh7s/HimbBkx096GNfE8i8AHwurtf1VhMubm5vmDBgv28IhGRr78bHm09SdjMizbFpJ3k5OQ8dz8qJo2JNCKeI2FfAJUEU4u1pQL1reFKA7KBP5vZnyNlCYCZ2W7gdHevO7VJZMTsX8BhMYtcREREJMbitibM3cuBPGBMnUNjCO6SrGstMAw4stbrQWBl5H1952DBnOURBAv+RURERFqkeI6EAfwWeNzM3gPeAi4D0gmSK8zsMQB3v9DdK4C6e4J9Duxy96W1ym4D3gVWAN2AqwiSsH3dcSkiIiISmrgmYe7+lJn1BH5OMN24lGBasThSpdH9whpwMPBHoBewmWAd2Sh3fy8GIYuIiIg0i3iPhOHufwD+0MCxE/dx7jRgWp2ya4FrYxOdiIiISHzo2ZEiIiIiIVASJiIiIhICJWEiIiIiIYgqCTOz75rZqbU+3xp5VNArZpYW+/BEREREWqdoR8KmVb8xs+HAzcDvgXbA3bELS0RERKR1i/buyL4Ez2oE+A7wN3efaWbzgVdiGpmIiIhIKxbtSNhOgodkA4wGXo2831yrXERERET2IdqRsDeBu81sEXAUMDFSfjjwWSwDExEREWnNoh0JuxIoJ0i+LnP3dZHycWg6UkRERKTJohoJc/c1wLfrKb8mZhGJiIiItAFR7xNmZh3NbKKZ3WhmB0fKBphZcuzDExEREWmdohoJM7OBBIvxDyJ4cPYzwJfA5ZHPl8Y6QBEREZHWKNqRsHuA+UAqsKNW+fPASbEKSkRERKS1i/buyJHAMe5eaWa1y1cD6TGLSkRERKSV259nR7arp6wPwV5hIiIiItIE0SZh84GptT67mXUDbgdeillUIiIiIq1ctNORU4HXzawQ6Ag8BQwENgDfjXFsIiIiIq1WtPuErTOzI4HvA8MJRtL+CPyPu+9o9GQRERERqRHtSBiRZOtPkZeIiIiI7Id9JmFmNh54wd0rIu8b5O5zYxaZiIiISCvWlJGwvwK9gM8j7xviQGIsghIRERFp7faZhLl7Qn3vRURERGT/RZVUmdkoM9srcTOzRDMbFbuwRERERFq3aEe2Xgfqe1D3wZFjIiIiItIE0SZhRrD2q66ewLYDD0dERESkbWjSFhVm9nzkrQNPmNmuWocTgaHA2zGOTURERKTVauo+YWWRnwZsAmpvzFoOLAJmxTAuERERkVatSUmYu18MYGargLvcXVOPIiIiIgcg2scW3d5cgYiIiIi0JU3ZMf9D4AR332RmS6h/YT4A7n5ELIMTERERaa2aMhI2B6heiN/YjvkiIiIi0kRN2TH/9vrei4iIiMj+02OIRERERELQlDVhja4Dq01rwkRERESapilrwrQOTERERCTGoloTJiIiIiKxoTVhIiIiIiHQPmEiIiIiIdA+YSIiIiIh0D5hIiIiIiGI6tmR1cxsAJAd+Vjg7p/ELiQRERGR1i+qJMzMegL/DZwFVP272F4ELnH3shjHJyIiItIqRXt35MPAQOB4oGPkNQroB8yKbWgiIiIirVe005GnAaPd/Z1aZW+Z2U+AV2MXloiIiEjrFu1IWCmwrZ7y7YCmIkVERESaKNok7D+Be8wso7og8v7uyDERERERaYL9eYB3P2CVma2NfM4AdgLfIFgzJiIiIiL7oAd4i4iIiIRAD/AWERERCYEe4C0iIiISgqiSMDNrb2a3m9nHZrbTzCprv5orSBEREZHWJtqRsF8CFxHcDVkF/D/gAYLtKa5oSgNmdoWZfRpJ4vLM7PhG6p5gZm+bWZmZ7TCzj8zs+nrqTTCz5Wa2K/LzO1Fel4iIiEhcRZuEfRe4zN0fAiqB59z9KuA2YMy+Tjaz7wH3AjOAXOBtYJ6Z9WnglK+A3xPsyj8Y+BVwu5nVJHxmdizwFPA/wJGRn8+Y2X9EeW0iIiIicRNtEpYKLI+8/wo4OPL+78CpTTh/KvCIu89y9wJ3nwKUAJfXV9nd89z9SXdf5u6fuvsTwCsEj02qdg3wurtPj7Q5Hfi/SLmIiIhIixRtErYaSI+8X0nwGCOAY4EdjZ1oZu2BEcD8OofmAyOb8uVmlhup+0at4mPrafOVprYpIiIiEoZonx35LDAaeJdgWvEvZjaJYMPW3+zj3EOARGBDnfINwCmNnWhma4CUSLy3u/uDtQ73aqDNXg20NRmYDJCWlsb7778PQHp6Op07d2blypUAdO/enf79+7N48WIAEhMTycnJobCwkG3bgic3ZWdns3HjRjZsCL6+d+/etG/fnqKiIgB69OhBnz59yM/PB6Bdu3YMGzaMgoICduwIctbBgwdTWlpKaWkpAH379sXMWLVqFQA9e/YkLS2NpUuXAtChQweGDBnCsmXL2LVrFwBDhw6lpKSEsrLgyVGZmZm4O8XFxQCkpKSQkpLC8uXBIGanTp3Izs5myZIlVFRUAJCTk8Pq1avZtGkTAP3796e8vJw1a9YAkJqaSnJyMgUFBQB06dKFrKws8vPzqawM7snIzc2lqKiIzZs3AzBw4EC2b9/OunXrqP7z7tatG4WFhQB07dqVww47jMWLF+PumBm5ubmsWLGCrVu3ApCVlcWWLVsoKSlRP6mf1E/N0E/Qg9airKwsJv0kEi/m7vuu1dDJZscQjDh97O4v7qNuOrAWOMHdF9YqvxU4392zGjm3H3AQcAxwJ3C1uz8eOVYOXOruj9WqfyEwy907NBZTbm6uL1iwYB9XKSLSet3waOtJwmZetCkm7SQnJ+e5+1ExaUykEdGOhO3B3d8lGBVrii8IFvOn1ilPBdbv43s+jbxdYmapwDTg8UjZ+v1pU0RERCRMUW/WambDzewxM/tX5PW4mQ3f13nuXg7ksfddlGMI7pJsqgSg9gjXOzFoU0RERCSuohoJM7PzgceABcDLkeJjgPfM7EeRuxcb81vgcTN7D3gLuIxgof+DkfYfA3D3CyOfpwCfAoWR80cB1wN/qNXmvcBCM7sJ+BvwHeAk4Lhork1EREQknqKdjpwO/MLdZ9QuNLOfEezh1WgS5u5PmVlP4OdAGrAUON3diyNV6u4XlkiwBiwT2A18AtxEJGmLtPm2mZ0X+f7/jNT5nrv/I8prExEREYmbaJOwFODpesqfAX7RlAbc/Q/sOZJV+9iJdT7fA9zThDb/Cvy1Kd8vIiIi0hJEuybsdeDEespPZM+9u0RERESkEfscCTOz8bU+zgPuMLOj+PddkccA4wnuWBQRERGRJmjKdGR903w1G57Wch8NTDOKiIiIyJ72mYS5e9TbWIiIiIhI45RgiYiIiIRgfzZrPcPMFprZF2ZWamZvmNnpzRGciIiISGsVVRJmZpcSPMT7E+BGgj27PgWeNbNLYh+eiIiISOsU7T5hNwJT3f3+WmX/bWZ5BAnZn2IWmYiIiEgrFu10ZB/g7/WUzwP6Hng4IiIiIm1DtEnYavZ+WDbAqUBxPeUiIiIiUo9opyPvAu4zs+HA25GybwEXAFNiGZiIiIhIaxZVEubuD5nZ58B1BLvkAxQA33X352IdnIiIiEhr1eQkzMySCKYdF7r7s80XkoiIiEjr1+Q1Ye6+G5gLdG2+cERERETahmgX5ucDA5sjEBEREZG2JNokbBpwt5mdY2aHmlly7VczxCcibdyrr77KN7/5TUaMGME999yz1/EHHniAY445huOOO45zzjmHzz77rObYxIkTyczM5LzzztvjnClTpnD88cdz3HHHcdFFF/HVV181+3WIiNQVbRL2EjCMYFpyFVAaeX0R+SkiEjOVlZXccMMNPP3007zzzjvMmTOHjz76aI86RxxxBAsWLGDRokWcddZZ3HbbbTXHpkyZwoMPPrhXu9OnT+fNN99k0aJF9O7dm4cffrjZr0VEpK5ot6g4qVmiEBGpR15eHv369SMzMxOA8ePHM2/ePAYNGlRT5/jjj695f9RRR/H000/XfD7hhBNYtGjRXu1269YNAHdn586dmFkzXYGISMOalISZWWfgN8A5QDvgVeAqd/+iGWMTkTaupKSEjIyMms/p6enk5eU1WP+JJ57glFNOaVLbP/3pT3n11VfJysril7/85QHHKiISraZOR94O/IhgOvIvBLvm/1czxSQiErWnn36axYsXM2VK0/aNfuCBB1i+fDmHH344zz6rXXdEJP6amoSNB37s7pPd/WrgDOAcM0tsvtBEpK1LS0tj7dq1NZ/XrVtHWlraXvX+7//+j7vvvpvZs2fToUOHJrefmJjI+PHjeeGFF2ISr4hINJqahB0KvFn9wd3fA3YD6c0RlIgIwPDhwykqKqK4uJjy8nLmzp3L2LFj96jz4YcfMnXqVGbPnk1KSso+23R3ioqKat7PmzePww47rFniFxFpTFMX5icC5XXKdkdxvohI1JKSkpg5cyYTJ06ksrKS888/n+zsbGbMmEFubi7jxo3jtttuY9u2bVx88cUA9O7dm9mzZwNw+umns2LFCrZt28aQIUP4/e9/z0knncQVV1zB1q1bcXeGDh3KXXfdFeZlikgbZe6+70pmVcD/ArtqFY8D3gC2Vxe4+1mxDrA55ebm+oIFC8IOQ0QkNDc82iPsEGJm5kWbYtJOcnJynrsfFZPGRBrR1OnIR4F1QFmt1xPAZ3XKRL629rUp6Ntvv82JJ55ISkoKzz235/Pqp02bxsiRIxk5ciRz586tKZ88eTLf/OY3GTlyJFdeeSUVFRXNfh0iIvL10KTpRHe/uLkDEQlT9aagc+fOJT09ndGjRzN27Ng99qPq3bs3DzzwAPfff5X4gQYAABEiSURBVP8e586fP5/8/HwWLlzIrl27OOusszjllFPo1q0b5557Lg899BAAkyZN4vHHH+eSSy6J67WJiEjLFO2O+SKtUu1NQdu3b1+zKWhtffr0YciQISQk7PnX5qOPPmLkyJEkJSXRpUsXBg8ezGuvvQbAmDFjMDPMjOHDh7Nu3bq4XZOIiLRsSsJEqH9T0JKSkiadO3ToUF577TW2b99OWVkZixYt2mNbBYCKigqefvppRo8eHdO4RUTk60t3N4ocoJNPPpnFixczduxYevbsydFHH01i4p5b6F1//fUce+yxHHvssSFFKSIiLY2SMBGaviloQ6677jquu+46IFj7NXDgwJpjd955J2VlZfzud7+LXcBfU7oTT0Tk3zQdKULTNgVtSGVlJRs3bgRg2bJlLFu2jJNOCp51/9hjj7FgwQJmzZq111oyERFp2zQSJkLTNgV9//33ueCCC9i8eTN///vf+fWvf80777xDRUUFp59+OgBdu3bloYceIikp+Kt13XXXceihh3LaaacBcOaZZ3LDDTeEdp0iItJyKAkTiRgzZgxjxozZo+zmm2+ueT98+HCWLVu213kdO3bk3XffrbfN0tLS2AYpIiKthuZHREREREKgJExEREQkBErCREREREKgJExEREQkBFqYL62S9qMSEZGWTiNhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiGIexJmZleY2admttPM8szs+EbqppnZbDP7yMwqzeyReur8yMy8nlfHZr0QERERkQMQ1yTMzL4H3AvMAHKBt4F5ZtangVM6AF8Avwb+0UjT24G02i933xmruEVERERiLd4jYVOBR9x9lrsXuPsUoAS4vL7K7r7K3a9y90eAjY206+6+vvYr9qGLiIiIxE7ckjAzaw+MAObXOTQfGHmAzXcys2IzW2NmL5pZ7gG2JyIiItKskuL4XYcAicCGOuUbgFMOoN1C4BIgH+gKXA28ZWY57r6ibmUzmwxMBkhLS+P9998HID09nc6dO7Ny5UoAunfvTv/+/Vm8eDEAiYmJ5OTkUFhYyLZt2wDIzs5m48aNbNgQXFLv3r1p3749RUVFAPTo0YM+ffqQn58PQLt27Rg2bBgFBQXs2LEDgMGDB1NaWkppaSkAffv2xcxYtWoVAD179iQtLY2lS5cC0KFDB4YMGcKyZcvYtWsXAEOHDqWkpISysjIAMjMzcXeKi4sBSElJISUlheXLlwPQqVMnsrOzWbJkCRUVFQDk5OSwevVqNm3aBED//v0pLy9nzZo1AKSmppKcnExBQQEAXbp0ISsri/z8fCorKwHIzc2lqKiIzZs3AzBw4EC2b9/OunXrqP7z7tatG4WFhQB07dqVww47jMWLF+PumBm5ubmsWLGCrVu3ApCVlcWWLVsoKSmJqp9ak9bUT61J9f87vu5/n6BHc/9RxU1ZWVlM/r8nEi/m7vH5IrN0YC1wgrsvrFV+K3C+u2ft4/wXgS/c/Uf7qJcIfAC87u5XNVY3NzfXFyxY0MQrkK+TGx5tPb9YZl60KewQYkb90vKoT/aWnJyc5+5HxaQxkUbEc03YF0AlkFqnPBWI2Roud68E/gUcFqs2RURERGItbkmYu5cDecCYOofGENwlGRNmZsARBAv+RURERFqkeK4JA/gt8LiZvQe8BVwGpAMPApjZYwDufmH1CWZ2ZORtN6Aq8rnc3ZdHjt8GvAusiNS5iiAJq/eOSxEREZGWIK5JmLs/ZWY9gZ8T7Oe1FDjd3YsjVerbL2xxnc/fBoqBzMjng4E/Ar2AzZH6o9z9vdhGLyIiIhI78R4Jw93/APyhgWMn1lNm+2jvWuDamAQnIiIiEid6dqSIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASZiIiIhICJSEiYiIiIRASVhIXn31Vb75zW8yYsQI7rnnnr2O79q1i0suuYQRI0ZwyimnsHr16j2Or1mzhkMPPZT77ruvpiwnJ4dvfetbjBo1ipNPPrnZr0FERET2X1LYAbRFlZWV3HDDDcydO5f09HRGjx7N2LFjGTRoUE2dJ554goMPPpi8vDzmzJnDtGnT+NOf/lRz/JZbbmH06NF7tf3888/Ts2fPuFyHiIiI7D+NhIUgLy+Pfv36kZmZSfv27Rk/fjzz5s3bo87LL7/MeeedB8DZZ5/NwoULcXcAXnrpJfr27btH0iYiIiJfL0rCQlBSUkJGRkbN5/T0dEpKShqsk5SURLdu3di4cSNfffUV9957LzfccMNe7ZoZEyZM4KSTTuKRRx5p1msQERGRA6PpyK+ZO++8k8svv5yDDjpor2Mvv/wy6enplJaWMn78eA4//HBGjhwZQpQiIiKyL0rCQpCWlsbatWtrPq9bt460tLR662RkZLB79262bNlCcnIyeXl5PP/880ybNo3NmzeTkJBAx44dmTRpEunp6QCkpKRwxhlnkJeXpyRMRESkhdJ0ZAiGDx9OUVERxcXFlJeXM3fuXMaOHbtHnXHjxvHkk08C8Nxzz3H88cdjZrz88svk5+eTn5/PZZddxrXXXsukSZPYtm0bW7duBWDbtm28/vrrZGdnx/3aREREpGk0EhaCpKQkZs6cycSJE6msrOT8888nOzubGTNmkJuby7hx4/jhD3/IZZddxogRI+jRowcPP/xwo22WlpZywQUXALB7924mTpzIKaecEo/LERERkf2gJCwkY8aMYcyYMXuU3XzzzTXvO3bsuM/F9TfddFPN+8zMTN58882YxigiIiLNR9ORIiIiIiHQSFgM3PBoj7BDiJmZF20KOwQREZE2QSNhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISAiVhIiIiIiFQEiYiIiISgrgnYWZ2hZl9amY7zSzPzI7fR/0TIvV2mlmRmV12oG2KiIiIhC2uSZiZfQ+4F5gB5AJvA/PMrE8D9fsBL0fq5QJ3APeZ2YT9bVNERESkJYj3SNhU4BF3n+XuBe4+BSgBLm+g/mXAOnefEqk/C3gUuP4A2hQREREJnbl7fL7IrD2wHfi+uz9Tq/wBYKi7n1DPOQuBJe7+01pl5wKzgc6A7Uebk4HJkY9ZQGEMLi8eDgG+CDsI2Yv6peVRn7RMX6d+6evuKWEHIa1fUhy/6xAgEdhQp3wDcEoD5/QCXq2nflKkPYu2TXf/I/DHJkfdQpjZv9z9qLDjkD2pX1oe9UnLpH4R2ZvujhQREREJQTxHwr4AKoHUOuWpwPoGzlnfQP3dkfZsP9oUERERCV3cRsLcvRzIA8bUOTSG4I7G+rzTQP1/uXvFfrb5dfW1m0JtI9QvLY/6pGVSv4jUEbeF+VCzncTjwBXAWwR3P/4YGOLuxWb2GIC7Xxip3w9YCswCHgK+BfyBYCH+nKa0GbeLExEREYlCPKcjcfenzKwn8HMgjSDBOr1WstSnTv1Pzex04HcEW06sA66qTsCa2KaIiIhIixPXkTARERERCejuSBEREZEQKAkTERERCYGSMBFpNczMav8UEWnJtCasBTOz3sBAgv3QqoBCd9f+ZyJNVJ2Muf5HJyItkJKwFsrMLgcuAXKAbcBKYA3wLvA3dy80swR3rwoxzDbFzDq5+46w45C9mVkCcDaQQvBc2bXAG+7+eaiBiYg0QklYCxTZcmMlcDfwXwS/WE4BTgSyCZKxa919uZmZ/pXf/MysB5APvAQ8Abxd/edeuw/MbBCwzt23hBZsG2NmXYH/Bk4iGDFeAziwE3gDeNzdP9Lflfgxs3ZAP6DY3XeFHY9IS6U1YS3TD4CP3f1X7l7m7h+5+/3uPhH4CcG/9F80s0P0SyVufkjwOKwRwEJgpZn9p5ll1UrADgX+QvBweYmfq4Asgv0BU4HzgXuAJcCpwEwzS9Hflbj6KbAYeNDMvm1mvcwssXYFM+tmZuMiCZtIm6QkrGUqB7qa2VAAM+tgZu0B3H0RwS+ZnQS/YCQ+jgD+DJwJ5AJPA98HlpvZu2Y2mSBRO8zdi8ILs00aCzzq7v8EiPyj5QngSuA6gtHjx0OMry36HvAewZrWvxE8gu43ZnacmXWP1PkBcJu7V4QUo0jolIS1TH8lmFa5xsy6uvsudy+PrHvB3VcDXwK9wwyyrTCzDsBy4DN3/9zdP3T3nwFHAadFjk0DpgN3hhZoG2RmSQRPyZhgZimRssTIeslKd19I8Ciz3maWE2asbUWkHyqAWe5+PNCXYLr4TIJR5AVmdiNwDfCP0AIVaQG0JqyFqXVr/dnAvUAywajLHwiG93sDowjWig1z91UhhNnmRBKxHu6+PjKt4rVvijCzE4EFQB93XxNSmG2SmR0D/A/BP15+6+4b6hw/FCgAstx9bQghtilmlgacByx391fqHMsFLo0c7wEcqj6RtkxJWAtlZgcTPEtzJPAdgoeXA6wn2LLicXefFk50bUv1gm4z6w9sq/1LvtaxW4EfuXv/8CJteyKjwwnAxcAMgufhzgGeAlYTTCOfCQx296PDirOtMbNOBP9Q2Vl7z7Za6yenE6zhyw0rRpGWQElYC2Jm3wAuIFjH8gWwg2DacRHB1hTtCNZY/N3dPw4rzrakVp9MBT4HdgMlwDPAXHffFvklM4ngrsgXQwu2jYv8w+VHBGuNjgS2Eqyd/Cdwh7tr6iuOGrob1cw6A+8Df3Z3Td9Lm6YkrAUxs0eAIcALwEaCqchhwOEECcDP9Yskvhrok1xgEMFWCL9x9/mhBdiGmVk3YGvtX/SRkbGOwEHAUIKRS/2diZP6+qSeOh0JFu7/xd3L4xacSAukJKyFiIymbCUYol9Yq6wP8B8E6yj6A9919/dDC7QNaaRPegPHEIx+9QW+rz6JPzN7iOAOvPcI9qPaa282M+vh7pu0R1h8NLFPDnb3L+MenEgLpLsjW47BwKcE21MAwfoJdy9296eBbxNMTZ4bUnxtUUN98pm7P0Ow1mgr6pO4M7PvEyTBdwPPEWx/MN7MBkbWI2FmBwF/NrNhSsCaXwN98h0zG1CrTzoBj1ZvvyPS1mkkrIWI/M/pRYKNWC8EPqn7SCIzmwL82N2PDCHENkd90nKZ2SygEpgJjAcuAgYAhcDLwGsEG7je6+7tw4qzLVGfiERPI2EtROSZhLcAnYDHgAvN7NDIv+arF7OeQLAnksSB+qRliuwN9inwpbsXuftd7j4MOJrgMUUXEWzrch/apDUu1Cci+0cjYS1MZJj+F8BZBA/ufgcoJXh2ZAlwqbsvCS/Ctkd90vJEnuWZGnkmZHugos4C/e8RPEJquLt/EFacbYn6RCR6SsJaqMjWCGcA5xDcZr8UeMbdPwo1sDZMfdKyRe6MNHevNLNJBNNencOOqy1Tn4g0TknY10DkESxV+64p8aI+adnMbCqQ6O6/CTsWCahPRPamJExEWh0zawdUKlFuOdQnIntTEiYiIiISAt0dKSIiIhICJWEiIiIiIVASJiIiIhICJWEiIiIiIVASJiIiIhICJWEiIiIiIfj/jjwg7z22NVYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_histogram([cal_results.get_counts('cal_11')], legend=[r'Preparing in $|11\\rangle$'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Correcting Using the Calibration Matrix\n",
    "\n",
    "If we have the calibration matrix $\\mathbf{A}$ which gives the transformation between the distributions, \n",
    "$$\\tilde{P}_{\\rho} = \\mathbf{A} \\cdot P_{\\rho}$$\n",
    "then to work back to $P_{\\rho}$ we just need to invert $\\mathbf{A}$,\n",
    "$$P_{\\rho} = \\mathbf{A}^{-1} \\cdot \\tilde{P}_{\\rho}$$\n",
    "However, there are some issues to watch out for.\n",
    "\n",
    "Below, we simulate a calibration with 200 shots and then prepare the state $|11\\rangle$ and measure with 2000 shots. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Corrected Results using Matrix Inversion:\n",
      "{'00': 91.36875000000059, '01': 46.12083333333413, '11': 2015.9520833333343, '10': -153.44166666666717}\n"
     ]
    }
   ],
   "source": [
    "# Generate the calibration circuits\n",
    "qr = qiskit.QuantumRegister(2)\n",
    "meas_calibs, state_labels = complete_meas_cal(qubit_list=[0,1], qr=qr)\n",
    "# Generate a noise model for the 2 qubits\n",
    "noise_model = noise.NoiseModel()\n",
    "for qi in range(2):\n",
    "    read_err = noise.errors.readout_error.ReadoutError([[0.9, 0.1],[0.25,0.75]])\n",
    "    noise_model.add_readout_error(read_err, [qi])\n",
    "backend = qiskit.Aer.get_backend('qasm_simulator')\n",
    "job_w_noise_cal = qiskit.execute(meas_calibs, backend=backend, noise_model=noise_model, shots=200)\n",
    "job_for_correction = qiskit.execute(meas_calibs[-1], backend=backend, noise_model=noise_model, shots=2000)\n",
    "cal_results = job_w_noise_cal.result()\n",
    "raw_results = job_for_correction.result()\n",
    "meas_fitter = CompleteMeasFitter(cal_results, state_labels)\n",
    "meas_filter = meas_fitter.filter\n",
    "print(\"Corrected Results using Matrix Inversion:\")\n",
    "print(meas_filter.apply(raw_results, method='pseudo_inverse').get_counts('cal_11'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you run the above code several times, you will likely see some negative counts. This is because there is some statistical noise in the calibration matrix which means that the corrected results will be unphysical. To correct for this we can find the $P_{\\rho}$ which is closest to reproducing the measured output,\n",
    "$$\\text{min}_{P_{\\rho}} ||\\tilde{P}_{\\rho}-\\mathbf{A}\\cdot P_{\\rho}||$$\n",
    "but where all elements of $P$ are non-zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Corrected Results using Constrained Least Square:\n",
      "{'00': 18.683235581467915, '01': 61.34919107554396, '11': 1919.9675733425333, '10': 4.5479099506544676e-10}\n"
     ]
    }
   ],
   "source": [
    "print(\"Corrected Results using Constrained Least Square:\")\n",
    "print(meas_filter.apply(raw_results).get_counts('cal_11'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Tensored Mitigation\n",
    "\n",
    "Tensored mitigation applies for local errors. In a fully-tensored mitigation all the errors act on a single qubit. In this case, the calibration matrix $\\mathbf{A}$ has the form\n",
    "$$ \\mathbf{A} = \\bigotimes_{i=1}^n A_i, $$\n",
    "where $A_i$ is a specific calibration matrix of the $i$th qubit.\n",
    "\n",
    "Fully-tensored mitigation requires only two calibration circuits, as opposed to $2^n$ circuits in the general case. Moreover, the calculation of the estimate of $P_{\\rho}$ does not require construction of the full calibration matrix $\\mathbf{A}$, thus saving the need to store an $2^n$ times $2^n$ matrix in memory; this is true for both pseudo-inverse and least-squares methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Noisy calibration matrices:\n",
      "[array([[0.91 , 0.228],\n",
      "       [0.09 , 0.772]]), array([[0.907, 0.253],\n",
      "       [0.093, 0.747]])]\n",
      "Corrected Results using Matrix Inversion:\n",
      "{'00': -0.41701417848210554, '01': 21.82374200722836, '11': 1916.592680280161, '10': 62.00059189109197}\n",
      "Corrected Results using Constrained Least Square:\n",
      "{'01': 21.607661075488608, '11': 1916.6337878419722, '10': 61.75855108253915}\n"
     ]
    }
   ],
   "source": [
    "# Generate the calibration circuits\n",
    "qr = qiskit.QuantumRegister(2)\n",
    "mit_pattern = [[0], [1]]\n",
    "meas_calibs, state_labels = tensored_meas_cal(mit_pattern=mit_pattern, qr=qr)\n",
    "# Generate a noise model for the 2 qubits\n",
    "noise_model = noise.NoiseModel()\n",
    "for qi in range(2):\n",
    "    read_err = noise.errors.readout_error.ReadoutError([[0.9, 0.1],[0.25,0.75]])\n",
    "    noise_model.add_readout_error(read_err, [qi])\n",
    "backend = qiskit.Aer.get_backend('qasm_simulator')\n",
    "job_w_noise = qiskit.execute(meas_calibs, backend=backend, noise_model=noise_model, shots=1000)\n",
    "cal_results = job_w_noise.result()\n",
    "meas_fitter = TensoredMeasFitter(cal_results, mit_pattern)\n",
    "print(\"Noisy calibration matrices:\")\n",
    "print(meas_fitter.cal_matrices)\n",
    "job_for_correction = qiskit.execute(meas_calibs[-1], backend=backend, noise_model=noise_model, shots=2000)\n",
    "raw_results = job_for_correction.result()\n",
    "meas_filter = meas_fitter.filter\n",
    "print(\"Corrected Results using Matrix Inversion:\")\n",
    "print(meas_filter.apply(raw_results, method='pseudo_inverse').get_counts('cal_11'))\n",
    "print(\"Corrected Results using Constrained Least Square:\")\n",
    "print(meas_filter.apply(raw_results).get_counts('cal_11'))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Tags",
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
